A comparative study of dynamical simulation methods for the dissociation of 

molecular Bose-Einstein condensates 
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We describe a pairing mean-field theory related to the Hartree-Fock-Bogoliubov approach, and ap- 
ply it to the dynamics of dissociation of a molecular Bose-Einstein condensate (EEC) into correlated 
bosonic atom pairs. We also perform the same simulation using two stochastic phase-space tech- 
niques for quantum dynamics — the positive P-representation method and the truncated Wigner 
method. By comparing the results of our calculations we are able to assess the relative strength 
of these theoretical techniques in describing molecular dissociation in one spatial dimension. An 
important aspect of our analysis is the inclusion of atom-atom interactions which can be problem- 
atic for the positive-P method. We find that the truncated Wigner method mostly agrees with the 
positive-P simulations, but can be simulated for significantly longer times. The pairing mean-field 
theory results diverge from the quantum dynamical methods after relatively short times. 

PACS numbers: 03.75.Nt, 03.65.Ud, 03.75.Gg, 03.75.Kk 



I. INTRODUCTION 

The dissociation of a molecular Bose-Einstein conden- 
sate (BEC) P, S i il into correlated atom pairs is a 
process analogous to parametric down-conversion in op- 
tics. Down-conversion involving photons has been piv- 
otal in the advancement of quantum optics by allowing 
for the generation of strongly entangled states. In the 
same way, molecular dissociation has emerged as an av- 
enue to generate strongly entangled ensembles of atoms 
in the field of quantum-atom optics. This matter-wave 
analog is of additional interest, however, as it gives rise to 
the possibility of performing tests of quantum mechanics 
with mesoscopic or macroscopic numbers of massive par- 
ticles rather than with massless photons. For example, 
the atom pairs formed during dissociation have Einstein- 
Podolsky-Rosen (EPR) type correlations in position and 
momentum, and one can envisage a demonstration of 
the EPR paradox with ensembles of correlated ultra-cold 
atoms 0, @, g. Also, molecular BECs can be formed 
by either two bosonic or two fermionic atoms; the lat- 
ter offers the possibility of a new paradigm in fermionic 
quantum atom optics. 

Experimental progress in the field of ultra-cold quan- 
tum gases has reached the stage where investigation of 
atom-atom correlations is now possible [H, Q . For exam- 
ple, in 2005 Greiner et al. Q measured atom-atom corre- 
lations resulting from the dissociation of molecules 
into fermionic atoms. Such advances have been achieved 
through the development of techniques for the measure- 
ment of noise in absorption images [3, S Hi fill and atom 
detection using microchannel plate detectors [12 |. 
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In this paper we consider correlations between bosonic 
atoms produced, for example, in the dissociation of ^^Rb2 
dimers. Whilst molecular dissociation of ^^Rb2 has been 
experimentally realised Q, atom-atom correlations have 
not yet been measured in these experiments due to the 
short molecular lifetimes. Experimental advances, how- 
ever, may soon result in the production of BECs of ro- 
vibrationally stable ground-state molecules [l^l , in which 
case the present analysis will become experimentally rel- 
evant. This paper serves to further the understanding 
of atom-atom correlations in the molecular dissociation 
process. Previous analytic and numeric work in this 
area has been restricted to the short time limit, where 
the effects of s-wave scattering interactions are negligible 
0, i, 0, [H, [H, [13, [H, [H, Ea, Jl, M ■ However , if a full 
quantitative description of atom-atom correlations is to 
be obtained, the effects of spatial inhomogeneity and s- 
wave scattering interactions on correlation strength must 
be addressed [a. Ell ■ To this end, we provide numerical 
results beyond the short time limit for the case of a spa- 
tially inhomogeneous molecular condensate with atom- 
atom interactions included in the model. 

The other contribution made in this paper is a com- 
parison of the performance of three simulation meth- 
ods describing the dynamics of BECs beyond Gross- 
Pitaevskii theory. Since the experimental realisation of 
Bose-Einstein condensation in 1995 [2^, the dynamics of 
weakly interacting BECs have often been successfully de- 
scribed by applying a mean-field theoretic approach lead- 
ing to the Gross-Pitaevskii equation (GPE) [U. How- 
ever, as the GPE neglects quantum fluctuations, its abil- 
ity to describe the full BEC dynamics is limited to cases 
where the effects of quantum fluctuations are negligible. 

Incorporating the effects of quantum fluctuations when 
modelling quantum many-body systems is necessary to 
describe, for example, the correlation dynamics which 
play a significant role in more recent experiments, such 
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as molecular dissociation. As a result of this, much ef- 
fort has been directed at developing theoretical methods 
that go beyond mean-field theory in their description of 
the dynamics of ultra-cold quantum gases [1^, |2^, \2^ . 
Several techniques have been used in the analytical and 
numerical investigation of the dissociation of a molecular 
BEC and the atom-atom pair correlations resulting from 
this process. For instance, dissociation can be treated an- 
alytically using the undepleted, classical molecular field 
approximation for the case of uniform condensates [31; 
a more recent development is the analytic treatment of 
nonuniform condensates using a perturbation theory in 
time [2I]. As the name suggests, the undepleted molec- 
ular field approximation assumes that the number of 
molecules remains constant throughout the dissociation 
process. Hence, it is only valid for short dissociation 
times when depletion is negligible, corresponding to a 
conversion of < 10% of the molecules into atoms [l4.[29|. 
Although useful in some circumstances, the obvious lim- 
itations of the analytic treatment beyond this regime ne- 
cessitates an alternative approach. 

In this paper we compare three simulation techniques 
using molecular dissociation as an example: two stochas- 
tic phase-space methods known as the positive-P (30l.[3]| 
and truncated Wigner [s^ . Issj methods, and a pairing 
mean-field theory known as the Hartree-Fock-Bogoliubov 
(HFB) method [13, [H, [H, [13. The positive-P rep- 
resentation method provides an exact quantum treat- 
ment of the dissociation problem for inhomogcneous sys- 
tems, with s-wave scattering interactions and molecu- 
lar depletion incorporated. Extensive work has been 
conducted usin g th e p ositive P-reprcsentation method 
P, [1, M M MM, m,MM, with savage et al. ^, ^ 
in particular, analysing both position and momentum 
pair-correlations in molecular dissociation. 

Unfortunately, the positive-P approach is also lim- 
ited to relatively short simulation times. For exam- 
ple, when one neglects the atom-atom interactions com- 
pletely, the positive-P simulations are successful only for 
durations corresponding to about 50% conversion [l^ . 
For typical experimental systems, divergent trajectories 
and large sampling errors arise during the evolution so 
that the problem becomes intractable beyond this time 
scale (40l.l4l|. The problem becomes worse when one in- 
cludes atom-atom s-wave scattering interactions; in this 
case the dissociation durations that can be simulated us- 
ing the positive-P method are limited to only ~ 5%, and 
at best 10%, conversion [3l. This prevents one from us- 
ing the positive P-representation method to determine 
the effects of s-wave scattering on the atom-atom corre- 
lation strength over time. Due to this limitation, there 
is a subsequent lack of knowledge regarding the effects 
of s-wave scattering on correlation dynamics for realistic 
condensates. This motivates further numerical investi- 
gation of atomic correlations in molecular dissociation 
using approximate methods, and to this end we consider 
the truncated Wigner and HFB methods to elucidate the 
relative performance of the methods in the context of 



molecular dissociation. 

This paper is structured as follows. Section II describes 
the system we have studied. Sections HI and IV provide 
an outline of the three simulation methods used in this 
work, including the relevant evolution equations and ap- 
proximations. Furthermore, it presents justification for 
the use of the three methods, discusses their inherent 
limitations and motivates the need for a comprehensive 
comparison of their relative performance. In Section V, 
we detail our work based on simulations of the coupled 
atom-molecule system, describing molecular dissociation 
in one dimension (ID). Finally, Section VI provides an 
overview of work extending the truncated Wigner simu- 
lations beyond short time scales. 

II. SYSTEM AND HAMILTONIAN 

We consider a molecular BEC which is dissociated into 
pair-correlated atoms by way of a magnetic Feshbach 
resonance. The quantum field theory effective Hamilto- 
nian describing this coupled atom-molecule system can 
be written as [3, [l3l , 



H = I d^l Yl *l(x)ffo,.(x)*.(x) 

i—a,m, 

+ E ^^!(x)^J(x)*,(x)*,(x) 

i.j—a.m 

+ ^(^L(x)«'^(x)+i/.c.)| (1) 

where ^a.m(x, t) are the atomic/molecular field oper- 
ators that annihilate an atom or molecule at position 
x. The field operators satisfy the commutation relation 
[*i(x, t), ^](x',t)] = (5iji5(x-x'). The atomic/molecular 

free-particle Hamiltonians, _ffo,a(x) and i/o.m(x), are 
given by, 

^0,a(x) = -^Vl + hVai^), (2) 

2nia 

i?o,m(x) = -^\7l + hV„,{x) + 2h\A\, (3) 

where TOq is the atomic mass and m„j = is the 

molecular mass. The atomic/molecular trapping poten- 
tials are given by Va(x) and Kn(x) = 214 (x). The de- 
tuning A in Eq. ([3]), or dissociation energy 2fi,|A|, corre- 
sponds to an overall energy mismatch of 2Ea — Em be- 
tween the free atom states 2Ea at the dissociation thresh- 
old and the bound molecular state Em- Hence, the pro- 
cess of dissociation begins with an initially stable, Em < 
2Ea, molecular BEC and a magnetic field sweep onto the 
atomic side of the Feshbach 

negative detuning A), resulting in the formation of atom 
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pairs. For a molecule at rest, the excess dissociation en- 
ergy is converted into the kinetic energy of atom pairs, 
which for the most part will posses s equal and opposite 
momenta ±ko, where fco = |ko| = ■\/2toq| A|/fi,. 

Returning to Eq. ([T]), Uij represents the two-body s- 
wave interaction strengths for atom-atom, atom-molecule 
and molecule-molecule scattering events. For example, 
Uaa = 'iTrhas/ma where is the atomic scattering length 
{qs = 5.4 nm for ®^Rb). The term x is the atom- molecule 
coupling and is responsible for coherent conversion of 
molecules into atom pairs, where the mechanism for con- 
version is via a Feshbach resonance [13, [3, . How- 
ever, in appropriately chosen rotating frames the equa- 
tions can easily be recast for conversion via optical Ra- 
man transitions [itI , [38| . In our numerical work the atom- 
molecule coupling remains switched on for the total evo- 
lution time. Also, we assume that the trapping potentials 
are switched off when the coupling x is switched on at 
t = and the evolution occurs in free space. 

The Hamiltonian ([1]) conserves the total number of 
atomic particles 

N = 2{N^{t)) + {Na{t)) = const, (4) 

with Ni{t) = /dx*J(x, i)*j(x, t) {i = a,m) and N, = 
(Ni). We begin our simulations with the molecular BEC 
in a coherent state and the atomic field in the vacuum 
state, and so iV = 2(iV„(0)). 



III. STOCHASTIC METHODS FOR BEC 
DYNAMICS 

After being developed in the field of quantum optics, 
phase-space methods have been successfully applied to 
matter- wave physics and have been used in many studies 
of the quantur n dy namics of complex many-body systems 

such as BECs [3, [H [H [11 [H, H, M M M M M M 

[52l . [53l . [S^ l . Phase-space representation methods rely on a 
mapping between the quantum operator equations of mo- 
tion and the Fokker-Planck equation (FPE) which in turn 
can be interpreted as a set of stochastic differential equa- 
tions (SDEs). Two distributions commonly used for this 
purpose can be traced back to the Glauber-Sudarshan P- 
distribution and the Wigner distribution [i^, [5^ . Along 
with the HFB method, phase-space techniques are cen- 
tral to this paper and hence will be discussed briefly in 
order to develop a context for the numerical results pre- 
sented in Sec. V and VI. 



A. The Positive P-representation Method 

The positive P-representation method [10, [3l|, [H, [10, 
|49| enables one to perform first-principles calculations of 
the quantum dynamics of multi-mode quantum many- 
body systems, including BECs. It relies on exploiting 
the positive P-represcntation of the density matrix, for 



which there exists a mapping between the master equa- 
tion and a set of c-number SDEs that can be solved nu- 
merically. The stochastic trajectory averages calculated 
using the positive P-representation method correspond 
to the normally-ordered expectation values of quantum 
mechanical operators [ij]. If stochastic sampling er- 
rors remain small during the time evolution, any observ- 
able can, in principle, be calculated using the positive-P 
method. 

The positive-P approach requires one to double 
the phase-space by defining two independent complex 
stochastic fields ^^(x, and $i(x, {i = a,m) corre- 
sponding to the operators ^i{x,t) and 5'|(x, t), respec- 
tively O] , with ^f* (x, t) ^ (x, t) except in the mean. 
Using the Hamiltonian in Eq. ([1]), the stochastic differ- 
ential equations describing the quantum dynamical evo- 
lution are, in the appropriate rotating frame [H, [l3] , 



dt 2ma 



+ A/-ix*mCl + V-*t^ma*a*m/2(C2 + <3) 



dt 2ma 



+ V^tW^A72(C6-<7) + V^C^^^^Cio- (5) 

Here the Cj{x, t) = 1, 2, . . . , 10) are real, independent, 
Gaussian noises with (^^(x, t)) = and correlations in 
time and space given by {Q{x,t)(h{^' ,t')) = SjkS{x — 
x')S{t-t'). 



B. The Truncated Wigner Method 

The truncated Wigner method is another useful phase- 
space technique for describing the quantum evolution of 
a Bosc- Einstein condensate [s^, [11] . Unlike the positive 
P-representation method it is an approximate method as 
it involves neglecting (or truncating) third-order deriva- 
tive terms in the evolution equation for the Wigner func- 
tion. This is necessary in order to obtain an equation 
in the form of a Fokker-Planck equation which can then 
be mapped onto a stochastic differential equation. The 
third-order terms can, in principle, be represented via 
stochastic difference equations, however, these are more 
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unstable than the positive-P equations [56| . The advan- 
tage of the truncated Wigner method hes in the inclusion 
of initial quantum noise, allowing the model to incorpo- 
rate quantum corrections to the classical field equations 
of motion and treat a different set of problems to a Gross- 
Pitaevskii equation or other classical field approaches. 
Although it has been shown that the truncated Wigner 
approach can give erroneous results, particularly for two- 
time correlation functions js^], it can be accurate for a 
wide range of problems pro vided the particle density ex- 
ceeds the mode density |46l [ssj . 

It can be shown that using the truncated Wigner ap- 
proximation (TWA) and the Hamiltonian in Eq. ^ , the 
stochastic difl[erential equations governing the dissocia- 
tion are, in the appropriate rotating frame, 



dt 



ih 
2m„ 



9*™ ^h .v2vi,„_,^[/_|vl/,|2,I,„^_,|vl,2. (6) 



dt 



2m,i 



Whilst these equations are deterministic, quantum fluc- 
tuations are included by way of a noise contribution in 
the initial state for the molecular and atomic fields. The 
addition of this initial vacuum noise ensures that the 
initial state of \E',„ and represent the Wigner func- 
tion of an initial coherent state BEG and an initial vac- 
uum state, respectively. The respective stochastic aver- 
ages with the Wigner distribution function correspond to 
symmetrically-ordered operator products, so that the cal- 
culation of observables represented by normally-ordered 
operator products needs appropriate symmetrisation. 



IV. PAIRING MEAN-FIELD THEORY FOR BEG 
DYNAMIGS 



in the model. However, they are neglected in our work 
as they are negligible on the time scales under considera- 
tion. This is one of the main differences between the HFB 
and truncated Wigner approaches, as the latter includes 
molecular fluctuations. By including the fluctuation op- 
erator, Xa(x), the atomic fleld is treated to higher order 
than the molecular field in our HFB formalism. This 
is necessary as the atomic fluctuations play an intrin- 
sic dynamical role in the molecular dissociation process 
and also allow one to consider atomic pair correlations. 
Finally, it is assumed that the initial molecular state is 
a coherent state and any expectation values of greater 
than two atomic fluctuation operators are factorised us- 
ing Wick's theorem [6^, thereby assuming that the quan- 
tum state of the system is Gaussian. 

With these approximations and our Hamiltonian for 
the system, given in Eq. we can derive a set of cou- 
pled PDEs for the atomic mean-field 0a (x), molecular 
mean- field </>,„ (x) and the first-order correlation functions 
Ga(x, x') and GAr(x, x'). Solving these coupled evolution 
equations one can then model the dynamics of dissocia- 
tion of a molecular BEG. 

Within the HFB formalism the evolution equations de- 
scribing the molecular dissociation process are given by 



^'/'a(x) 

dt 



dt 



ih 
2m„ 



-iUaa [|(/.a(x)|2-t-2Gw(x,x)] M^) 
-iUaaGA{yi,^)4>*a{yi) - iX0m (x)(?!); (x) , 



(7) 



ih 
Arua 



V2(/,,„(x) -2i|A|0,„(x) 



-if/mm|0m(x)p(?ir„(x) - i^[0^(x) 

+Ga(x,x)], 



(8) 



Pairing mean-field theory - as a simplified version of 
the Hartree-Fock-Bogoliubov (HFB) approach ^59|, iG^, 
[6l| - has been applied to the problem of molecular dis- 
sociation in Refs. [l^, [2^, although these works only 
considered spatially uniform systems. Our present HFB 
study extends the analysis to nonuniform condensates 
and represents the third method we use in describing 
the dynamics of the molecular dissociation. This ap- 
proach involves an approximation to the full quantum 
evolution retaining only the lowest order atomic fluctua- 
tions. More precisely, one writes the atomic fleld operator 
^'a(x) in terms of the atomic mean- field (/)a(x) = (^a(x)) 
and the lowest order atomic fluctuations ^^(x), such 
that, *I'a(x) = '/'a(x) + Xa(x). Thc atomic fluctua- 
tions can be approximately represented by their lowest 
order correlation functions, the normal and anomalous 
densities, GAr(x, x') = {x\i^Xaip^)) ^^^^ G/i(x, x') = 
(Xa(x')Xa(x)), respectively [lllel. 

In our implementation the molecular fleld is treated 
as a mean-field, with 0m (x) = (5'm(x)). As suggested 
in Refs. [29, molecular fiuctuations can be included 



aGA(x,xO 
dt 



--{mx')x{x),H]) 



9GAr(x, x') 
di 



|0a(x)|^ + |0„(x')r + GAr(x,x) 



+ GAr(x',x') Ga(x,x') 



0a(x)2G^(x,x') +0a(x')'Gjv(x,x') 



+ Ga (x, x) G^ (x, x' ) + Ga (x' , x') Gjv (x, x' ) 



ia(x)2+GA(x,x) <5(x-x') 



0™(x)[G^(x,x')+<5(x-x')] 



+ 0m(x')G7v(x,x') 



(9) 
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-j^V^G^.(x,x')-^V^,G^.(x,x') 



2iUa 



|0a(x)|^-|0a(x')r+Gw(x,x) 



Gw(x',x') Gw(x,x') 



- iUaa 



0,(x)2g:i(x,x') + 0:(x')'Ga(x,x') 



Ga(x, x)G:i(x, x') - G*A^', x')Ga(x, x') 



,(x)G:i(x,x')-C(x')Ga(x,x') 



(10) 



where GAr(x, x) is the density of the noncondensed 
atoms. This follows from the expression for the total 
density of atoms, ('i'a*a) = |0a|^ + (xlxa), where \(j)a\'^ 
is the density of the condensate atoms and {x\xa) is the 
density of the noncondensed atoms. 

In our simulations we neglect (jja (x) as the atomic field 
does not develop since all the terms on the right-hand 
side of Eq. ^ are multiplied by the field or its conjugate. 
The physics here is similar to that of a non-degenerate 
optical parametric oscillator with phase diffusion [6^ . 
More particularly, only the sum of the phases of each 
correlated atom pair is known, fixed to the phase of the 
molecular BEC, whilst the relative phase is unknown and 
takes an arbitrary value. It follows that the individual 
phases of these correlated modes are also arbitrary and 
consequently, no common phase to the atomic field exists 
across the entire range of momenta. 

The potential role of the HFB and truncated Wigner 
methods, despite being approximate techniques, is to 
model realistic inhomogeneous condensates in which the 
effects of s-wave scattering interactions on atom-atom 
pair correlations can be quantified and compared with 
experimental data. Moreover, both methods present the 
possibility of describing physics in regimes for which the 
positive-P representation method is computationally in- 
tractable. Prior to this work, however, there has been no 
attempt to undertake a comprehensive comparison of the 
performance of all three of these methods when applied 
to the same problem, although there has been compar- 
isons of the positive-P and truncated Wigner methods 
when investigating BEC collisions [6^. This motivates 
the comparative study of these approximate methods and 
the positive-P approach. 

There are further potential advantages in developing 
the HFB method for application to such problems. The 
positive-P representation and truncated Wigner methods 
require the averaging of many trajectories (corresponding 
to quantum mechanical ensemble averaging) and there- 
fore requires multiple runs. In contrast, since the HFB 
method is not a stochastic technique it only necessitates 
a single run but at the expense of the dimensions of the 
problem being doubled. Also, in many ways it is a more 
intuitive method, with the derivation highlighted here 
being an extension of the well-known Gross-Pitaevskii 
approach. 



V. COMPARISON OF POSITIVE-P, HFB AND 
TRUNCATED WIGNER RESULTS 



A. Parameter values 

We now present ID simulations for the dissociation of 
a ^'^Rb2 molecular BEC with rua = 1.44 x 10~^'^ kg and 
nim = ^rUa- For computational simplicity we consider 
an effective one-dimensional (ID) system by assuming 
strong harmonic confinement in the transverse direction. 
All parameters are chosen to be close to a typical exper- 
imental system, with the exception of a relatively small 
value for the detuning |A| so that the computational grid 
need not be too large. In practice, the detuning should 
be such that the total dissociation energy 2h\A\ is much 
larger than the thermal energy due to the finite temper- 
ature of the system; here, we assume a zero-temperature 
condensate. At the same time the detuning |A| should 
be smaller than the frequency of the transverse harmonic 
trap potential, so that transverse excitations are sup- 
pressed and the dynamics of dissociation remains in ID. 

We set Umm = and Uam = in our simulations; it 
is the role of atom-atom s-wave scattering that is of par- 
ticular interest in this work, and to this end we perform 
simulations with both Uaa = and Uaa 7^ 0. Setting 
Umm to be zero is unlikely to be entirely physical, how- 
ever for a more realistic value we find that the positive-P 
simulations become intractable after a very short time, 
making our goal of a comparison impossible. Addition- 
ally, we find that there is no significant difference between 
the TWA simulation for Umm — and Umm 7^ 0, and so 
this has no practical implications for our study. 

Similar considerations apply to the atom-molecule in- 
teractions which are set to Uam =0. At the mean-field 
level, the atom-molecule interactions in the equations of 
motion for the atomic field would initially appear as an 
effective spatially varying detuning that depends on the 
molecular BEC density profile; these interactions can be 
neglected if the total dissociation energy 2h\A\ is much 
larger than the respective interaction energy per atom. 
For our choice of A and the molecular BEC peak density 
(see below) this would in turn require an atom-molecule 
scattering length of < 0.1 nm. For more realistic values 
of the atom-molecule scattering length (assumed to be in 
the few nanometers range) the approximation would re- 
quire absolute detunings in the kHz range or higher and 
it would improve with increasing |A|. 

The initial molecular BEC density is taken to be gaus- 
sian, 



,(x,i = 0) =noe-"'/'^' 



(11) 



corresponding to a trapping frequency of 0.15 Hz in the 
x-direction with a harmonic oscillator length of 50 /xm, 
where ng = 1.83 x 10^ m~^ is the molecular BEC peak ID 
(linear) density. The size of the one-dimensional quanti- 
sation box was chosen to be P = 6.5 x 10"'' m and the 
lattice grid was composed of 512 points. The atom-atom 
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interaction strength is given by Uid = Uaa/'^A = 2iu±as 
[1^, where A = nl^ is the confinement area in the trans- 
verse direction, with l± = ^/h/miJ± being the trans- 
verse ground-state harmonic oscillator length and uj± is 
the transverse oscillation frequency. The atom-molecule 
couphng in ID is given by xiD — x-\/27rZj[ = 6.7 x 10^'^ 
m^/^ s~^ [1^, and is switched on for the total evolu- 
tion ti me. The detuning A = —258 s^^ and hence, 
ko = y/2ma\A\/h ~ 8.41 x 10^ m^^ is the resonant mo- 
mentum at ±fco. 

With the reduction of the coupling constants to their 
ID counterparts, the equations of motion in previous sec- 
tions are unchanged, except that all propagating fields 
(and the respective noise sources in the positive-P equa- 
tions) are now understood as ID fields, while the operator 
Vx is replaced by d"^ /dx^. 

In all our simulations, we assume that the atomic field 
is initially in a vacuum state and that the molecular con- 
densate is initially in a coherent state, with density profile 
given by Eq. (fTTjl . The trapping potential is turned off 
when the dissociation coupling xiD is switched on, with a 
Feshbach sweep into the dissociation regime A < 0. Dis- 
tinct from the implementation of the stochastic methods, 
the HFB simulations assume that the molecular conden- 
sate remains in a coherent state during the dynamical 
evolution. Also, the initial atomic fluctuation flelds are 
assumed to be Ga(x, x') = GAr(x, x') = and the molec- 
ular fluctuations are omitted. 

The numerical codes for solving the evolution equa- 
tions for the methods under consideration, were imple- 
mented using the XMDS simulation package [63|. AH 
stochastic simulations were performed for the case of 
10,000 trajectories. In the following sections, we have 
verified the accuracy of the results presented by ensur- 
ing, for instance, invariance of results for different lattice 
size and time step. Furthermore, we were able to perform 
benchmarking with the analytic result within the unde- 
pleted, molecular field approximation up until t = 0.06s 
(~ 10% molecular conversion), and more importantly, 
with the exact positive-P results. 



B. Initial comparisons 

We first perform simulations neglecting atom-atom 
interactions with Uid — for an initial number of 
molecules Nm{Q) = 1-62 x 10'^. We observe the formation 
of peaks in the atomic density at momenta ±fco as the 
dissociation energy (excess potential energy) is converted 
into the kinetic energy of the correlated atom pairs, with 
equal but opposite momentum ±fco- We verify that the 
value of fco agrees with the predicted value, given in 
SecES 

In Fig. [T] we provide a plot of the total fractional num- 
ber of molecules Njn{t)/Njn{0) and the total fractional 
number of atoms Na{t)/2Nm{0) as a function of time t, 
normalised to the total molecule A'„i(0) and total atom 
number 2Nm{0), respectively. Although this result does 



not include the effects of s-wave scattering, it does al- 
low one to compare the performance of the methods. It 
can be seen that all the methods agree until t ~ 0.14 
s, which corresponds to the conversion of ^ 10% of the 
molecules. It is found that whilst the truncated Wigner 
method does extremely well when compared to the ex- 
act results provided by the positive-P method, the HFB 
method diverges substantially at longer times. 



1.2 



g 0.6 

E 

^ 0.4 



0.2 



■'■ molecules (HFB) 
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—molecules (+P) 
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t(s) 
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FIG. 1: (Color online) Comparison of the fractional particle 
numbers for t final ~ 0.20 s, for the case of a non-uniform con- 
densate with UiD = 0, for positive-P (red solid and dashed 
lines), truncated Wigner (black ■ and ►) and HFB (blue 
dashed lines) methods. The fractional atom numbers (ini- 
tially lower curves) and the fractional molecule numbers are 
shown. In this figure and throughout this paper, the positive- 
P and truncated Wigner results are for the case of 10,000 
trajectories. The error bars are shown and are essentially the 
thickness of the data lines in all figures. In all simulations per- 
formed the initial number of molecules is Nm{0) = 1.62 x 10^. 

The positive-P method becomes intractable at < 0.20 
s, and hence the ability to compare all three methods 
ceases beyond this point. Looking forward, when one in- 
corporates s-wave scattering the positive-P method will 
fail sooner [l3| and hence the truncated Wigner and HFB 
methods may be able to access a regime otherwise inac- 
cessible to numerical simulations for realistic non-uniform 
condensates. 



C. Observation of Phase Diffusion Processes 
During Dissociation of a Molecular BEC 

In this section we consider non-uniform condensates 
with s-wavc scattering interactions included. In Fig. [2l 
we plot the fractional particle number throughout the 
evolution, for the same parameters as in Fig. [1] but with 
scattering included. We choose an interaction strength 
of UiD = go = 1.04 X 10~^uj±as, which corresponds to 
^^Rb with transverse confinement of LLj±/2n = 30 Hz and 
s-wave scattering length of = 5.4 nm. From these re- 
sults it can be seen that the positive-P method fails be- 
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yond approximately iinax(+^') = 0.18 s, whilst the trun- 
cated Wigner method produces results beyond tmax(+^') 
and still does well in comparison to the positivc-P re- 
sults up to <max(+P)- As cxpectcd we find that 
the positive-P method fails for even shorter times as 
the interaction strength is increased. For example, with 
an interaction strength of Um = 32(7o we find that the 
positive-P method fails for fmax(+-P) ~ 0.05 s. 
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t (+P) 



0» " »■>■> " > " >■ " ►■'»• '*' 

1) 0.05 



0.1 
t(s) 



0.15 



0.2 



FIG. 2: (Color online) As in Fig. [T] except for the case of 
a non- uniform condensate with Uid ~ go- The positive-P 
method becomes intractable beyond tmax{+P) ~ 0.18 s, com- 
pared with f ~ 0.2 s when s -wave scattering is neglected. We 
also find that the number of molecules converted into atom 
pairs is decreased when s-wave scattering interactions are in- 
cluded and attribute this to phase diffusion. 



—ko — ko — > — 2fcoi Eind are observed in the positive-P and 
truncated Wigner results. They do not arise in the HFB 
results as the method does not allow for uncondensed 
molecules outside the initial condensate mode. 



D. Analysis of Atomic Pair-Correlation Functions 

We have also investigated atomic pair-correlations re- 
sulting from the dissociation process, for realistic non- 
uniform condensates including the effects of s-wave scat- 
tering. The strength of atom-atom correlations can be 
quantified using Glauber's second-order correlation func- 
tion 5(2)(k,k',t) [zl. 



5(^)(k,k',i) 



(a+(k,t)at(k',t)a(k',t)a(k, t)) 
(n(k,t))(fi(k',t)) ■ 



(13) 



with the momentum-space density at k given by n(k, t) = 
(n(k, i)) = (a^(k, i)d(k, i)) where the momentum-space 
field amplitudes are represented by the lattice-discretized 
momentum components a^(k) and a(k), which corre- 
spond to the continuous Fourier transforms of the fields 
in the limit Afc [l3|. This pair-correlation function 
describes the ratio of the probability of the joint detec- 
tion of atom pairs with k and k' to the product of the 
probabilities of independent atom detection at k and k'. 
From this it follows that g'''^\k,'k' ,t) = 1 for uncorrc- 
lated atom pairs, g^^^k, k', t) = 2 for thermally bunched 
atoms and ^'^^^(k, k', t) > 2 for strongly correlated atoms 



Here we also observe the signature of phase diffusion 
during molecular dissociation [6^, [t^I ■ By consid- 

ering Eq. ([5]) we are able to estimate the characteristic 
diffusion time for the dissociation process. 



2C/iz5(«'t(2; = 0,irf)*(.T = 0,id))' 

and verify that the process of phase diffusion is responsi- 
ble for suppressing molecular conversion, and hence, de- 
creasing the number of atom pairs formed. We also per- 
formed simulations for increased values of the atom-atom 
interaction strength, Uid = 2^0 and Uid ~ 32go. The 
results show that molecule conversion decreases with in- 
creasing interaction strength and further supported the 
order of magnitude estimates of the diffusion time. Un- 
fortunately, as in Sec. IV B[ we see that the HFB method 
still fails to adequately describe the dynamics of the 
molecular dissociation process for longer times, with the 
particle numbers only providing bounds for the true val- 
ues. Another feature indicative of the limitations of 
the HFB method is its inability to predict the forma- 
tion of peaks in the molecular momentum spectrum at 
±2fco |7l| . in addition to the main resonant momenta 
peaks formed at ifep- These secondary peaks arise due 
to atom-atom recombination processes fco + fco ^ 2fco and 




It . . . i 

0.05 0.1 0.15 0.2 

t(s) 

FIG. 3: (Color online) Plot of the atomic pair-correlation 
functions for back-to-back and collinear scattering processes, 

(2) (2) 

denoted gggiko, —ko,t) and (?^^(fco, fco, t). Results are shown 
for t final = 0.20 s, for a non-uniform condensate with Uid = 
0. The positive-P results (red solid and dashed lines), the 
truncated Wigner results (black ■ and ►) and the HFB re- 
sults (blue o and *) are shown. The collinear (red dashed, 
black ►■ and blue *) and the back-to-back pair-correlations 
(red solid, black ■ and blue o) are shown. 

We quantify pair-correlations arising due to momen- 
tum conservation which are present between atoms with 
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equal but opposite momenta, and pair-correlations aris- 
ing due to quantum statistical effects [i.e. the Hanbury- 
Brown and Twiss (HBT) bunching] between atoms 
scattered in the same direction. The atomic pair- 
correlations function for back-to-back (BE) and coUinear 
(CL) scattering processes, are denoted ggg{kQ,—kQ,t) 

and (7p^(fco, fco, t), respectively. These quantities are 
shown in Fig. [3] for the case of no s-wave scatter- 
ing and in Fig. 0] with scattering incorporated. The 
collinear correlation indicates HBT thermal bunching 
with g^|(fco, fco) t) = 2 until t ~ 0.10 s in both cases. The 
back-to-back correlation is super-bunched due to strong 
correlations between atom pairs with equal but opposite 

(21 

momenta, with ggg{ko, — fco,i) > 2 for short times. Be- 
yond t 0.10 s we observe both the collinear and back- 
to-back correlations are approximately equal, fall below 
two, and approach the uncorrelated or coherent level of 
5^^^ = 1, for the truncated Wigner and positive-P re- 
sults. The HFB results, on the other hand, fail to predict 
where g'^' approaches the coherent state level as stimu- 
lated processes become important. Toward the end of 
the simulation, the back-to-back correlation drops below 

(2) (2) 

the collinear correlation, gggik^, —ko, t) < gj^j^{ko, ko,t). 
This effect becomes more severe with increasing values of 
UiD and is also noticeable in the HFB results. 




0.05 0.1 0.15 0.2 

t(s) 

FIG. 4: (Color online) As in Fig. [S] except for the case of 
a non- uniform condensate with Uid ~ go- The positive-P 
method becomes intractable beyond tmax{+P)^ 0.18 s, com- 
pared with t ~ 0.2 s when s -wave scattering is neglected. 
By comparison with Fig. O we see that the back-to-back and 
collinear pair-correlation strength degrades at an increased 
rate when s-wave scattering is included. 

From Figs. [3] and 2] it is again clear that the truncated 
Wigner method is most successful in describing molecu- 
lar dissociation, with the positivc-P method intractable 
at longer times. It should be noted that the truncated 
Wigner results for this correlation function are not shown 
prior tot 0.05 s, where sampling issues arise due to the 
small number of atoms per mode. However, once the sig- 
nal is significant the results agree with positive-P. As 
seen in Sec. V B and C, the HFB method fails to fully 
describe the dynamics for longer times as the molecular 



field deviates from the assumed coherent state. 



VI. SIMULATIONS BEYOND tmax(+P) 

The numerical results we have presented indicate that 
the HFB method is unsuitable for quantitative correla- 
tion studies of molecular dissociation beyond the regime 
of the positive-P simulations. The HFB method be- 
comes invalid as it assumes a mean-field coherent state for 
molecules for the entire simulation time [l^, [2^ . How- 
ever, once molecular depletion reaches ~ 80%, this as- 
sumption is no longer valid and the method becomes in- 
creasingly inadequate as the regime of complete depletion 
is reached. This assertion is further supported in Fig. [5l 
which provides a surface plot of the molecular density 
in position space beyond tmax(+P)- Here we begin to 
observe the development of ripple effects which coincide 
with the reduction of the molecular condensate density 
and it is unlikely that the approximation of the molecular 
field as a coherent state is still valid. In this regime the ef- 
fects of quantum fluctuations become increasingly impor- 
tant and hence, we cannot rely on the HFB method. To 
remedy this, the inclusion of molecular fluctuations Xm hi 
the HFB formalism could be one avenue for future work. 
It should be stressed that there is value in using the HFB 
method, as it lies between the crude undepleted, molecu- 
lar field approximation and an exact quantum treatment. 
For instance, the HFB ap pro ach is suited to high energy, 
sparsely occupied modes j73| . In such cases, fluctuation 
effects are largely insignificant and the HFB method is 
valuable. 




FIG. 5: (Color online) Molecular density in position space 
nm{x) [in units of m~^] as a function of time for the HFB 
results. The molecular peak density is given by no = 1-83 x 
lO'^m-V 

With the HFB approach found to be invalid beyond the 
realm of the positive-P simulations, we look at extending 
the simulations using the truncated Wigner approach. 
Fig. [5] and [7] repeat the analysis given in Sec. |V] but with 
the extension to t = 0.40 s. In Fig. [5] we again look at 
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the fractional particle numbers for the cases neglecting 
and including atom-atom interactions. Beyond t ~ 0.2 s, 
for the UiD — case we observe the effects of atom-atom 
recombination. This is apparent due to the slight increase 
in the number of molecules, and corresponding decrease 
in the number of atoms, until t ~ 0.3 s. With the effects 
of s-wavc scattering included, we again witness a phase 
diffusion process which is responsible for decreasing the 
rate of molecule conversion. 



1= 

\ 0.8 
^EO.6 
\ 0.4 
=-«0.2 




mol. U_1 D = g_0 
• atoms U_1 D = g_0 
■■■atoms U_1D = 
— mol. U 1D = 



0.1 



0.2 
t(s) 



0.3 
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FIG. 6: (Color online) Fractional particle numbers for tfi„ai ~ 
0.40 s, for the case of a non-uniform condensate with Uid ~ 
(black solid line and red dashed line) and Um = go (green 
crosses and blue dotted), for the truncated Wigner results. 
The fractional atom number (blue dotted and red dashed) 
and the fractional molecule number (black solid and green 
crosses) are shown. For the Um = case for t > 0.20 s, we 
observe an increase in the molecular population from atom- 
atom recombination. With the inclusion of s-wave scattering 
interactions we again see the effects of phase diffusion. 



In Fig. [7] we provide the truncated Wigner results 
for the back-to-back and coUinear pair-correlation func- 
tions for the cases where s-wavc scattering interactions 
arc neglected and included. In both cases, the back-to- 
back correlation is exceeded by the coUinear correlation, 

(2) (2) 

i.e. gg^{ko, —ko,t) < j^{ko , ko , t) , with the effect be- 
coming more dramatic with time. It is interesting to 
note that when atom-atom interactions are neglected the 
back-to-back pair correlation eventually turns into anti- 

(2) 

correlation, ie. gssi^O: ~ko) < 1, for sufhciently long 
times when the molecular depletion is large. As the 
atomic density increases, we see atom- atom recombina- 
tion which is not correlated at the two momenta consid- 
ered, so that atoms are not removed equally from each 
of the modes under consideration. Overall, by using the 
truncated Wigner method to go beyond the realm of the 
positive-P simulations, we are able to observe the effects 
of s-wave scattering on correlation dynamics for realistic 
inhomogeneous condensates. 
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FIG. 7: (Color online) Plot of the atomic pair-correlation 
functions for back-to-back (black solid and green crosses) and 
coUinear (blue dotted and red dashed) scattering processes, 
g^gg{ko, —ko,t) and g'^l^{ko, ko,t), respectively. Results are 
shown for t final = 0.40 s, for a non- uniform condensate with 
UiD — (blue dotted and green crosses) and Um — go (black 
solid and red dashed). The back-to-back correlation drops 
below the coUinear correlation. 



VII. CONCLUSIONS 

In this work we have compared three different theoret- 
ical approaches to the problem of dissociation of molec- 
ular Bose-Einstein condensates. We have considered the 
case where the atoms resulting from this dissociation pro- 
cess are not trapped, but move away from the parent 
molecules with momenta that are a function of the detun- 
ing. In particular, we have calculated atomic and molec- 
ular populations and analysed the effects of atom-atom 
interactions beyond the short time limit for inhomoge- 
neous condensates. We have also investigated quantum 
correlations, providing quantitative results for the back- 
to-back and coUinear pair-correlations, which cannot be 
calculated in the standard mean-field Gross-Pitaevskii 
approach. This is a subject of immediate interest as ex- 
periments which can measure these correlations can now 
be performed, particularly with metastable helium Q. 

In principle, the preferred theoretical method would be 
the stochastic integration of equations in the positive-P 
representation, as these give complete access to all prop- 
erties of the interacting many-body quantum system. In 
practice, however, the problems inherent in the integra- 
tion of these equations, especially when s-wave interac- 
tions of any appreciable strength are present, mean that 
the positivc-P equations are only useful for short times. 
Another method which has been widely applied to model 
BEG dynamics is the HFB approach. In some sense this 
is equal to the commonly used linearisation procedures of 
quantum optics, and similarly to that area, we find that 
we must be careful with its validity. In fact, we have 
shown here that the HFB approach will sometimes be- 
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come inaccurate on shorter time scales than those which 
give problems in the positive-P representation approach. 
Although it does present computational advantages in 
that the equations need only be solved once by contrast 
with the phase-space representations where averages need 
to be taken over many realisations, we see that it is also 
not useful for all parameter regimes. This could be reme- 
died, at least in part, by including the effects of molecular 
fluctuations. However, this is a cumbersome and compu- 
tationally expensive process. 

We have found that the most useful of the methods 
is the truncated Wigner representation. Although the 
approximations necessary to obtain stochastic differen- 
tial equations mean that the mapping from the quan- 
tum Hamiltonian is not exact, we find that the truncated 
Wigner method agrees with the first-principle positive-P 
results whenever such a comparison is possible to make, 
ft also has the advantages of not suffering from the sta- 



bility problems of the positive-P representation and is 
valid over longer times than the HFB approach. In con- 
clusion therefore, we find that while the positive-P and 
HFB approaches are useful in some regimes, the trun- 
cated Wigner representation is best suited to this prob- 
lem. 
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